%==========================================================================
% LMJ_SHELL.m 
%
% This controls loading for most programs to be run 
% It also provides an initial call to the model 
% Following variables must be in memory 
%
% cucd                   % Current Working directory 
% onscreen=1;            % Will query to continue
% flag_makelocalcopy =1; % Copy the settings and loading file to the outpath 
% root ='orig';          % Folder 
% spec ='log';           %    SubFolder 
% subf ='april 2010 pr1';%        Subfolder 
% set_fname=[];          % Name of settings file, default is settings 
% add2outpath=[];        % Append to LOADPATH name before the current date
% 
% Subfunctions used: 
% LOAD_LMJFILES 
% LMJ_MODFEVAL 
%
% Created      : April 6 2010
% Last modified: 
% June 2010   Allow for Missing data provided flag_mixfreq and flag_gtv 
%            defined 
% 
% TO DO 12/27/2011
% Remove flag_GTV 

% =========================================================================
cucd=pwd; 

% -------------------------------------------------------------------------
% 1. ROOTPATH from where will load settings file 
% -------------------------------------------------------------------------
loadpath= cr_dirpath(cucd,root,spec,subf); 
cd(loadpath); 
if isempty(set_fname); 
    settings; 
    set_fname='settings'; 
else 
    eval(set_fname); 
end 
cd(cucd); 
disp('Loaded settings'); 

%%=========================================================================
% Check settings for mixed frequency 
% i)    FLAG_MIXFREQ must equal 0 or 1 
% ii)   FLAG_GTV     must equal 0 or 1 
% =========================================================================
disp('Checking variable FLAG_MIXFREQ (0,1) exists'); 
if exist('flag_mixfreq','var')==0 
    error('FLAG_MIXFREQ, binary (0,1) not defined') 
else 
    if flag_mixfreq~=0 && flag_mixfreq~=1 
        error('FLAG_MIXFREQ, should only be (0,1)') 
    end 
    switch flag_mixfreq
        case 0 
            disp('FLAG_MIXFREQ==0, regular data, single frequency') 
        case 1 
            disp('FLAG_MIXFREQ==1, mixed frequency, data can have nans') 
    end            
    if flag_mixfreq==1
        disp(' '); 
        disp('Checking variable FLAG_GTV (0,1) exists');
        if exist('flag_gtv','var')==0
            error('FLAG_GTV, binary (0,1) not defined')
        end
        if flag_gtv~=0 && flag_gtv~=1
            error('FLAG_MGTV, should only be (0,1)')
        end
        switch flag_gtv
            case 0
                disp('FLAG_GTV==0, G time invariant')
            case 1
                disp('FLAG_GTV==1, G time varying, cummulator defined')
        end
        addsol.flag_gtv=flag_gtv; 
    end
end 
addsol.flag_mixfreq=flag_mixfreq; 

%% Check if this is a structural model or not. 
% If not, several outputs will not be produced. 
if exist('flag_structural','var')==0 || isempty(flag_structural) 
    flag_structural=1; 
    disp('Structural Model'); pause(0.2); 
else 
    disp('This is NOT a structural model'); 
    disp('Limited Output will be produced'); 
    pause(0.5); 
end 



%% -------------------------------------------------------------------------
% 2. Create model handle 
% -------------------------------------------------------------------------
funcmod=str2func(['mod',root,spec]); 

% -------------------------------------------------------------------------
% 3. Load priors, starting value, etc. 
% -------------------------------------------------------------------------
load_lmjfiles; 

if ~isempty(data_fname);
%     if sum( any(isnan(Y)) )
%         if flag_mixfreq==0
%             error('Y contains NANs but FLAG_MIXFREQ==0')
%         end
%     else
%         if flag_mixfreq==1
%             error('Y does not contain NANs but FLAG_MIXFREQ==1')
%         end
%     end
else 
    if onscreen==1
    disp('Not loaded any data'); 
    quer('c'); 
    end 
end 
if flag_valload==1
    dispaj('Starting value or mode loaded in PARLOADVAL'); 
end
if flag_prparload==1
    dispaj('Prior for parameters loaded'); 
end
if flag_prssload==1
    dispaj('Prior for steady state loaded'); 
end
if flag_gridload==1
    dispaj('Lower and Upper Grid for parameters loaded in PARGRID'); 
end
if flag_strvalsload==1;
    dispaj('Matrix of parameter starting values loaded in PARSTRVALS'); 
end

% To be completed 
% if flag_dataload==1;
%     dispaj('Data matrix DATAMAT loaded')
% end

% -------------------------------------------------------------------------
% 4. Create OUTPATH and make local copies of the file loaded and the
% settings file 
% -------------------------------------------------------------------------
%%
if isempty(add2outpath)
    outpath=cr_dir(loadpath,strdate); 
else 
    outpath=cr_dir(loadpath,[add2outpath,' ',strdate]); 
end 
dispaj('Created OUTPATH:',outpath); 

if flag_makelocalcopy==1 
    cd(loadpath); 
    [flag_success,message     ] =copyfile([load_fname,'.xls'],outpath); 
    if flag_success==1; 
        dispaj('Copied ',load_fname,' to OUTPATH'); 
    else 
        disp(message);         
        error('Cannot copy file loaded to OUTPATH'); 
    end     
    %temp_newname=[load_fname,' ',subf(2:end),' ',strdate,'.xls'];
    %temp_newname(find(isspace(temp_newname)==1))='_'; 
    %eval( ['!rename ',load_fname,'.xls ',temp_newname] );
    %load_fname=temp_newname; 
    %clear temp_newname; 
    
    
    [flag_success,message     ] =copyfile([set_fname,'.m'],outpath); 
    if flag_success==1; 
        dispaj('Copied ',set_fname,' to OUTPATH'); 
    else
        disp(message);         
        error('Cannot copy SET_FNAME to OUTPATH'); 
    end     
    clear flag_sucess message; 
            
end 
cd(cucd); 

% -------------------------------------------------------------------------
% 5. Check Options for NL Solver exist, if not assign default values
% -------------------------------------------------------------------------
disp(' '); 
if exist('solveopt','var')==0
    disp('Structure Options for NL solver not defined, assign default values');
    quer('c'); 
    solveopt=optimset('MaxFunEvals',200000,'Display','off',...
        'MaxIter',20000,'TolFun',1e-12);
else
    disp('Structure Options for NL solver exists');
end
dispaj(sprintf('TolFun      = %3.10g',solveopt.TolFun));
dispaj('MaxFunEvals =',solveopt.MaxFunEvals  );
dispaj('MaxIter     =',solveopt.MaxIter      );
dispaj('Display     =',solveopt.Display     );
if ~isempty( solveopt.TolX);
dispaj(sprintf('TolX        = %3.10g',solveopt.TolX));
end
if onscreen==1 
    quer('c'); 
end 

% -------------------------------------------------------------------------
% 6. Initial call to the function handle with either: 
% 1) Loaded starting value PARLOADVAL
% 2) First  in the matrix of starting values PARSTRVALS 
% Obtain SSNAMES 
% -------------------------------------------------------------------------
if flag_valload==1    
    disp('Solving model with loaded parameter value'); 
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames,ny,nx,nz,stateNames,shockNamesStru]=lmj_modfeval(funcmod,parloadval,solveopt,addsol); 
    numpar=size(parloadval,1); 
elseif flag_strvalsload==1; 
    disp('Solving model with first of starting valuess');
    [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames,ny,nx,nz,stateNames,shockNamesStru]=lmj_modfeval(funcmod,parstrvals(:,1),solveopt,addsol);  
    numpar=size(parstrvals(:,1),1); 
else
    error('Did not load starting value'); 
end 
if ~isequal(eu,[1;1]); 
    error('No Solution'); 
end 
if onscreen==1 
    quer('c'); 
end 

% =========================================================================
% Extracting position of SS parameters with prior
if flag_prssload~=0;
    display('extracting position of SS parameters with prior')
    ssposest=cellposition(deblank(prssnames),ssnames);
    printcell([ssnames(ssposest) num2cprec(ssvec(ssposest))]);
else
    display('Not loading PRIOR for SS');
    ssposest=[];
end

% =========================================================================
%% Check that if will change handles, models are identical and states exist
%% 

if exist('add2func','var')==1 && ~isempty(add2func) 
    dispaj('Switching function handle from ',func2str(funcmod)); 
    dispaj('to  ',[func2str(funcmod),add2func]);
    disp(' ');
    disp('Verifying solutions match...');
    pause(0.25);
    if flag_valload==1
        [~,~,difvar,difc]=checkmodels_lmj(funcmod,...
            str2func([func2str(funcmod),add2func]),parloadval,solveopt,addsol);        
    else
        [junk,junk,difvar,difc]=checkmodels_lmj(funcmod,...
            str2func([func2str(funcmod),add2func]),parstrvals(:,1),solveopt,addsol);
    end
    if difvar < 1e-5;
        disp('Solutions identical')
    else
        disp('Solutions do not match')
    end    
    if difc < 1e-5;
        disp('Constants identical');
    else
        disp('Constants do not match');
    end
    clear var1 var2 difvar difc
    if flag_valload==1
        [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames,ny,nx,nz,stateNames,shockNamesStru]=...
            lmj_modfeval(str2func([func2str(funcmod),add2func]),parloadval,solveopt,addsol);
    else 
        [GG, RR, CONS, eu, SDX, ZZ, initss, ssvec, flag, ssnames,ny,nx,nz,stateNames,shockNamesStru]=...
            lmj_modfeval(str2func([func2str(funcmod),add2func]),parstrvals(:,1),solveopt,addsol);
    end
end 
shonames=shockNamesStru.short; 

%% Check variables to report state_sel, spec_rep , shocks_sel exist 
% disp('Checking state_sel is correct....'); 
% junk=cellposition(state_sel(:),stateNames(:)); 
% if any(junk==0);error('State_Sel is incorrect');end 

disp('Checking specrep is correct....'); 
% junk=cellposition(specrep,stateNames(:)); 
% if any(junk==0);error('specrep is incorrect');end 
% if exist('shocks_sel','var')==1 && ~isempty(shocks_sel)
%     disp('Checking shocks_sel is correct...'); 
%     junk=cellposition(shocks_sel,stateNames(:)); 
% end 
if size(shonames)~=size(SDX,1); 
    error('Mistmatch Size SDX and Shock Names') 
end 
if  flag_mixfreq==1;     
    ZZ  =[ZZ;initss.AA];
end 
znames=cell(size(ZZ,1),1);
pos_obsinstates=zeros(length(znames),1); 
for ii=1:length(znames)
    temp      =find(ZZ(ii,:,1)~=0 ) ;
    if length(temp) > 1
        error('Need 1 entry only for each row of ZZ')
    end
     znames(ii)=stateNames( temp );
     pos_obsinstates(ii)=temp; 
end 
disp('Observables'); 
printcell([znames(:) Ynames(:)]); 
%quer('c'); 
clear junk znames initss; 


disp(' '); 
disp('End LMJ_SHELL'); 
disp('________________________'); 